% TEST_STOKES_ANNULUS_GRADUEX: gradient of the exact solution.

function gu = test_stokes_annulus_graduex (x, y)
uxx = @(x, y) (2.*(x-y).*y.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
	       (4+5.*x.^4-8.*x.^3.*y+4.*x.*y.*(5-4.*y.^2)+ ...
		5.*y.^2.*(-3+y.^2)+3.*x.^2.*(-5+6.*y.^2))+ ...
	       4.*x.*(x-y).*y.*(-4+x.^2+y.^2).* ...
	       (x.^5-8.*y-2.*x.^4.*y+20.*y.^3-6.*y.^5+ ...
		2.*x.^2.*y.*(5-4.*y.^2)+x.^3.*(-5+6.*y.^2)+ ...
		x.*(4+5.*y.^2.*(-3+y.^2)))+ ...
	       4.*x.*(x-y).*y.*(-1+x.^2+y.^2).* ... 
	       (x.^5-8.*y-2.*x.^4.*y+20.*y.^3- ...
		6.*y.^5+2.*x.^2.*y.*(5-4.*y.^2)+x.^3.*(-5+6.*y.^2)+ ...
		x.*(4+5.*y.^2.*(-3+y.^2)))+ ...
	       2.*y.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
	       (x.^5-8.*y-2.*x.^4.*y+20.*y.^3-6.*y.^5+2.*x.^2.*y.* ...
		(5-4.*y.^2)+x.^3.*(-5+6.*y.^2)+ ...
		x.*(4+5.*y.^2.*(-3+y.^2))));

uxy = @(x, y) (2.*(x-y).*y.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
	       (-8-2.*x.^4+12.*x.^3.*y+60.*y.^2-16.*x.^2.*y.^2- ...
		30.*y.^4+2.*x.^2.*(5-4.*y.^2)+x.*(10.*y.^3+10.*y.*(-3+y.^2)))+...
	       4.*(x-y).*y.^2.*(-4+x.^2+y.^2) .* ...
	       (x.^5-8.*y-2.*x.^4.*y+20.*y.^3-6.*y.^5+2.*x.^2.*y.* ...
		(5-4.*y.^2)+x.^3.*(-5+6.*y.^2)+ ...
		x.*(4+5.*y.^2.*(-3+y.^2)))+ ...
	       4.*(x-y).*y.^2.*(-1+x.^2+y.^2).* ...
	       (x.^5-8.*y-2.*x.^4.*y+20.*y.^3-6.*y.^5+2.*x.^2.*y.* ...
		(5-4.*y.^2)+x.^3.*(-5+6.*y.^2)+ ...
		x.*(4+5.*y.^2.*(-3+y.^2)))+ ...
	       2.*(x-y).*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
	       (x.^5-8.*y-2.*x.^4.*y+20.*y.^3-6.*y.^5+ ...
		2.*x.^2.*y.*(5-4.*y.^2)+x.^3.*(-5+6.*y.^2)+ ...
		x.*(4+5.*y.^2.*(-3+y.^2)))- ...
	       2.*y.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
	       (x.^5-8.*y-2.*x.^4.*y+20.*y.^3-6.*y.^5+ ...
		2.*x.^2.*y.*(5-4.*y.^2)+x.^3.*(-5+6.*y.^2)+ ...
		x.*(4+5.*y.^2.*(-3+y.^2))));

uyx = @(x, y) (-2.*(x-y).*y.^2.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
               (10.*x.^3+10.*x.*(-3+x.^2)-8.*x.^2.*y+2.*(5-2.*x.^2).*y+12.*x.*y.^2-4.*y.^3)- ...
               4.*x.*(x-y).*y.^2.*(-4+x.^2+y.^2).* ...
               (4+5.*x.^2.*(-3+x.^2)+2.*x.*(5-2.*x.^2).*y+(-5+6.*x.^2).*y.^2-4.*x.*y.^3+y.^4)- ...
               4.*x.*(x-y).*y.^2.*(-1+x.^2+y.^2).* ...
               (4+5.*x.^2.*(-3+x.^2)+2.*x.*(5-2.*x.^2).*y+(-5+6.*x.^2).*y.^2-4.*x.*y.^3+y.^4)- ...
               2.*y.^2.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
               (4+5.*x.^2.*(-3+x.^2)+2.*x.*(5-2.*x.^2).*y+(-5+6.*x.^2).*y.^2-4.*x.*y.^3+y.^4));

uyy = @(x, y) (-2.*(x-y).*y.^2.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).*...
               (2.*x.*(5-2.*x.^2)+2.*(-5+6.*x.^2).*y-12.*x.*y.^2+4.*y.^3)-...
               4.*(x-y).*y.^3.*(-4+x.^2+y.^2).*...
               (4+5.*x.^2.*(-3+x.^2)+2.*x.*(5-2.*x.^2).*y+...
                (-5+6.*x.^2).*y.^2-4.*x.*y.^3+y.^4)-...
               4.*(x-y).*y.^3.*(-1+x.^2+y.^2).*...
               (4+5.*x.^2.*(-3+x.^2)+2.*x.*(5-2.*x.^2).*y+...
                (-5+6.*x.^2).*y.^2-4.*x.*y.^3+y.^4)-...
               4.*(x-y).*y.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).*...
               (4+5.*x.^2.*(-3+x.^2)+2.*x.*(5-2.*x.^2).*y+(-5+6.*x.^2).*y.^2-...
                4.*x.*y.^3+y.^4)+...
               2.*y.^2.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).*...
               (4+5.*x.^2.*(-3+x.^2)+2.*x.*(5-2.*x.^2).*y+(-5+6.*x.^2).*y.^2-4.*x.*y.^3+y.^4));

gu = zeros (2, 2, size(x,1), size(x,2));
gu(1, 1, :, :) = reshape (uxx (x,y), [1, 1, size(x)]);
gu(2, 1, :, :) = reshape (uyx (x,y), [1, 1, size(x)]);
gu(1, 2, :, :) = reshape (uxy (x,y), [1, 1, size(x)]);
gu(2, 2, :, :) = reshape (uyy (x,y), [1, 1, size(x)]);

end
